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Abstract 

We propose a novel path sampling method based on the Onsager-Machlup (OM) action by 
generalizing the multiscale enhanced sampling (MSES) technique suggested by Moritsugu and 
coworkers (J. Chem. Phys. 133, 224105 (2010)). The basic idea of this method is that the system we 
want to study (for example, some molecular system described by molecular mechanics) is coupled to 
a coarse-grained (CG) system, which can move more quickly and computed more efficiently than the 
original system. We simulate this combined system (original + CG system) using (underdamped) 
Langevin dynamics where different heat baths are coupled to the two systems. When the coupling is 
strong enough, the original system is guided by the CG system, and able to sample the configuration 
and path space more efficiency. We need to correct the bias caused by the coupling, however, by 
employing the Hamiltonian replica exchange where we prepare many path replica with different 
coupling strengths. As a result, an unbiased path ensemble for the original system can be found 
in the weakest coupling path ensemble. This strategy is easily implemented because a weight for 
a path calculated by the OM action is formally the same as the Boltzmann weight if we properly 
define the path "Hamiltonian" . We apply this method to a model polymer with Asakura-Oosawa 
interaction, and compare the results with the conventional transition path sampling method. 
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I. INTRODUCTION 

Molecular dynamics (MD) simulation [l|, |2| is a rigorous and versatile approach for in- 
vestigating dynamic trajectory of molecular systems by numerical integration of Newton's 
equation. Its applicability has been quickly expanded in time and length scales owing to 
the development of hardwares (general purpose ones^such as PC clusters, GPGPU, super- 
computers, or special 



n 

force fields (Amber [4| 
numerical algorithms 



purpose ones such as Anton [3[), convenient softwares with rehab le 
Gromacs [sl, CHARMM [ej, NAMD [71 and others), and efficient 

2|. However, because the MD step size should be close to the vibra- 
tional time scales of the molecular system (~ Ifs), conventional MD simulations have not 
been applicable to processes that are extremely rare, slow, and/or diffusive processes. While 
the MD simulation is able to explore fast structural fluctuations around a basin of potential 
(or free) energy surface, it often fails to sample interbasin hopping, which describes various 
phenomena such as protein folding and chemical reactions 8|, M\ ■ 

To overcome this problem, several methods have been proposed and tested over the last 
two decades (for reviews, see 10|, lll|). One of the pioneering works is transition path 
sampling (TPS) developed by Chandler, Dellago, Bolhuis and coworkers [sHlO]. Since then, 
there has been effort to improve the efficiency of TPS calculations, such as transition interface 
sampling (TIS) [12| andpartial TIS 13[. Related to these methods are milestoning \\J^ and 
forward flux sampling 15|]. Studies on further improvement of the efficient path sampling is 
still under way. 



An alternative approach to path samp 



Onsager-Machlup (OM) action method 



search problems by Elber and coworkers 



16 



i ng is the action-based methods such as the 



17|. The OM action has been used for path 



18|, Eastman and coworkers [19|, and Orland and 
coworkers [20|]. In the early works, this method has been used to find the most dominant 
pathway, but recent advances have made it possible to explore the ensemble of paths 21 |. 
For instance, an efficient path sampling has been proposed previously 22| with the use of 
temperature replica exchange 23N25| . The computational advantage of the action-based 
methods is that the time step for the action integral can be taken to be as large as 100 fs (or 
more), which is much larger than a typical MD step (~ Ifs). Moreover, it is straightforward 
to carry out such calculations in parallel with respect to each time step and each path. 
This will certainly match the advantage of recent computational environment with hyper- 



parallel architectures such as GPGPU and supercomputer. Because of these computational 
demands, the ensemble of paths has not been easily obtained, but the fluctuations of paths 
as well as configurations should be important for biomolecules, which is worth pursuing to 
understand biological functions. 

In this paper, we suggest to enhance the efficiency of the OM action-based path sampling 



using the multiscale enhanced sampling (MSES) devised by Moritsugu and coworkers 26| 



The MSES method has been successfully applied to the calculation of free energy profiles 
of a small peptide (chigno 



in) [26|, intrinsically disordered protein (sortase) ^], and a 
protein-protein complex 28|]. In the original MSES, the system of interest is the canonical 
ensemble for the atomistic degrees of freedom, which will be called "MM DoF" . We assume 
a coupling between the MM DoF and a reference coarse-grained DoF, which will be called 
"CG DoF" , mimicking the typical behavior of the MM DoF in an economical manner. The 
total Hamiltonian of the combined system is defined as 

H = Hum + Hqg + A-^mmcg (1) 

where Hum, Hqg and Ai^MMCG represent the Hamiltonians of the MM and CG DoF, and 
the interaction between the MM and CG DoF, respectively. 

Assuming that the CG DoF is able to move around the CG configuration space quickly, 
the MM DoF will be dragged accordingly for a finite value of A, which leads to an efficient 
search of the MM configurations. However, our goal is to obtain the canonical ensemble 
with respect to MM-DoF, which is proportional to exp(— /Jifj^j^), in the absence of the MM- 
CG interaction, i.e. A = 0. To meet both of these ends, we use the Hamiltonian replica 



exchange 29|] among the combined MM-CG systems with different A values. In other words, 
we simulate a replicated set of the combined MM-CG systems with A^ > Ag ■ ■ ■ Xn^^ , where 
N is the number of replicas, and they are exchanged from time to time according to the 
Metropolis criterion [2j]. We then collect the MM ensemble with the smallest \n^^ (ideally 
^N^^ = 0), which should be equal to the unbiased canonical MM ensemble. 

Using the OM action-based methods, we can do the same as above: We consider an 
extended system where MM and CG DoF are coupled, and the OM action for the coupled 
system is defined as 

•S" = 5'mm + 'S'cG + ASmmcg (2) 

where Sum, Scg are the OM actions for the MM and CG DoF, and ASmmcg the interaction 



between them. Using A as the exchange parameter, the "canonical weight" for a MM path 
oc exp(— /J^mm) can be retrieved with the same logic as above. Note that in the formulation 
of MSES there is a subtle issue concerning the nonequilibrium nature of the combined 
system due to the coupling to different temperatures J26|, but this problem is circumvented 
by using the idea of "scaling" (see Appendix C). As shown below, when we only use the 
"conventional" OM method for the MM DoF (without any coupling to the CG DoF), we 
need a longer computation time for statistical convergence. Hence it is indispensable to use 
this type of method to enhance path sampling, especially when dealing with large molecular 
systems. 

This paper is organized as follows. In Sec. Ull we introduce the Onsager-Machlup action 
principle and formulate the multiscale enhanced sampling for path space using the OM 
action (MSES-OM). We also mention the model polymer system devised by Micheletti and 
coworkers 30|, |3l|, which is employed in the following calculations. In Sec. IIIIl we investigate 
the numerical performance of the MSES-OM method using the model polymer. We also 
examine the numerical accuracy of the method by comparing with the results using the 
conventional transition path sampling. In Sec. llVt we summarize the paper and discuss the 
further aspects of the MSES-OM method. 

II. FORMULATION 

A. Onsager-Machlup action for multi-dimensional systems 

Let us consider the case where the dynamics of the MM system is described by an 
overdamped Langevin equation, 

1 dU 



Xa = j^ + V^D^Vait) (3) 

where x^ is mass-weighted coordinates, U a potential energy, 7^ an intrinsic friction coeffi- 
cient. Da an diffusion coefficient, and rjait) a Gaussian- white noise satisfying {ria(t)ria'(t')) = 
6aa'S{t — t'). As Eq. ([3]) is a stochastic differential equation, the dynamics is characterized by 
the probability ("weight") of generating a path (or trajectory) x{t). According to Onsager 



and Machlup [32|, l33|, the weight is proportional to exp{—(3S[x{t)]) where ^[^(t)], which is 



called the Onsager-Machlup (OM) action, is a functional of x(t). In numerical calculations. 



the OM action is approximated by its discretized form as 

■'Vbcad — 1 



S[x{t)] 



AU 



+ E 



i=l 



■ M 9 



y , -^iXa,i+l - Xa,i) + Ves{{Xa,i}) 



a=l 



(4) 



where iVbead is the number of discretized segments of a path, or "beads" , M the number of 
the MM DoF (usually 3x number of atoms), AU = t/(xfin) — f^(a;ini) the potential difference 
between the initial and final states, and the effective "potential" is defined as 

^^ r 1 i.„T 1 

(5) 



Kff({xJ) = Yl 



a=l 






with 



Wn 



7q 



2At ■ ^^^ 

This is an effective frequency determined by the friction and the time interval for the path 
discretization AtoM- Here we used the notations U^ and U^x to represent the first and second 
derivatives of U with respect to x. 

To sample path space using the above OM action, we have employed the position- Verlet 
, l2| type scheme to integrate the Langevin equations, i.e.. 



Xn 



P'a, 



X„ 



Pa,i Ar 



m,; 



(1 - 7Ar/2K,, - 



dS 



dXr 



At + ^/2^kBTmArRaM^ 



'-'a, I — •*'q:,i 



(7) 
(8) 

(9) 



where Pa,i, mi, 7 are the fictitious momenta, mass, and friction coefficient, which are similar 
to those for the imaginary-time path integral simulations [SJI- Ri{t) is the Gaussian random 
variable with zero average and unit variance. Note that the time step to integrate the 
above equations of motion At has nothing to do with the time step for the discretization of 
the OM action AtoM- The numerical integration of these equations result in the canonical 
distribution exp{—pS). However, if the path space is huge, which is usually the case for 
complex systems, an efficient algorithm to sample the path space is required to reduce the 
computational load. This is the reason why we will introduce the MSES method 261-128 1. 



B. Multiscale enhanced sampling for path space 



In the formulation of MSES 26|, Moritsugu et al. introduced several replicas with dif- 
ferent coupling parameters Ai, A2, ■ • • , and exchange m-th and n-th parameters during the 
simulation of the combined system according to the probability 

Prnn = min(l, exp Amn) (10) 

where 

Amn = /3(Am — A„)(ifMMCG,m " -f^MMCG,n) (H) 

and P is the inverse temperature of the MM DoF. This is an application of the Hamiltonian 
replica exchange method 



29[ to the combined MM-CG system. The MM canonical ensemble 
is obtained from that of the combined MM-CG system with the smallest A. 

This formalism can be immediately applied to the problem of path sampling using the 
OM action. In this case the exchange probability is determined by 

Amn = /3(Am — A„)(S'MMCG,m — 'S'MMCG,n) (12) 

where the interaction OM action can be chosen as 

^MMCG = - 5^ miXo.,}) - xfy (13) 

i=l 

where 9i represents collective variables calculated from MM DoF {xa,i}, xf^ is the corre- 
sponding CG variables, and N^^ is the number of beads to discretize the CG path. The 
schematic picture of our MSES-OM method is shown in Fig. [H The basic strategy here is 
that the CG variables (blue, upper) move rather freely in path space to "guide" the MM 
variables (red, lower). 

Here it is important to notice on the scaling of the number of replicas required as the 



system size increases. As emphasized by Moritsugu et al. [261], the interaction "Hamiltonian" 
scales with the number of CG DoF not with that of MM DoF because it contains only CG 
DoF and the projected MM DoF onto the CG space. Hence this method does not suffer 
from the problems caused by the number of MM DoF as in the case of temperature replica 
exchange 2J, |25|. Furthermore, in this study, the number of CG beads is smaller than that 



of MM beads, i.e., iVbead (=240) is ten times smaller than A^bead (=2400), reducing the cost 
of calculations. 



In order for MSES to work efficiently, the CG DoF must move quickly and dominate 
the dynamics of the combined system. This situation is achieved, for example, when the 
total mass of the CG system is comparable to or larger than the MM system and the 
temperature of the CG system is much higher than that of the MM system. We thus take 
the CG mass = 20mo and CG temperature = SO/cbT, whereas mo is the "bead" mass of the 
model polymer (see Appendix A) and the temperature for the MM DoF is ksT c^ 300K. 
(Note that this CG mass is not a real mass because this is just introduced to carry out the 
"artificial" OM dynamics for the CG DoF.) Because two different temperatures are set for 
the MM and CG DoF, the resulting ensemble seems to exhibit a strongly nonequilibrium 
behavior. However, the desired information can be retrieved from the ensemble with the 
smallest coupling parameter X^^^ . The subtle point to use different temperatures for the 
MM and CG DoF are summarized in Appendix C. 

C. Model polymer 



In this section, we briefly explain the model system we use for path sampling, which is 
a coarse-grained polymer studied by Micheletti and coworkers 30|]. The potential energy of 
this model, U1 + U2 + f/3, consists of the Asakura-Oosawa (AO) type interaction f/3 and the 
finitely extensible elastic chain Ui + U2, 



U^ = e^^e-"("^^-"^^). 



where f, 



2r + r°. 



Uo 



U, 



l^iji dij 



i<j 

i 

(pksT 



1 ^r^ 



IQr^ 



\Ri — R 



J2 (2-~^. + 3r., - ^) 



f^e(f 



«J7' 



(14) 

(15) 
(16) 



■Jh ' ij 



Ri + Rj, 9(x) is the Heaviside function (=1 



for X > and for a; < 0), and the friction coefficient is given by 

7i = 67rr/i?i(l + 2.50). 



(17) 



The parameter values we employ here are almost the same as in the original paper 30 1 
except that we made a parameter four times larger, that is, we used a softer core. 

Micheletti et al. used the AO interaction to model the crowding effect due to RNA and 
proteins in a cell, and studied the (un) looping kinetics of the model polymer by changing the 
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FIG. 1: Schematic picture of the proposed multiscale enhanced samphng for path space using the 
OM action (MSES-OM). The blue (red) beads represent CG (MM) DoF, and the arrows between 
them with different thickness represent different interaction strengths. The thick arrows between 
the combined MM-CG systems schematically represent "exchange" due to Hamiltonian replica 
exchange (actually the coupling parameters are exchanged). 

size of the system. They found that the mean-first passage time (MFPT) for the unlooping 
dynamics does not depend on the size of the system, whereas the MFPT for the looping 

beads (five) , and study the 



3Q|. 



dynamics does. In this paper, we employ the smallest number o: 
unlooping dynamics, which takes place with ~ 3 /is time scale 

For path sampling using the OM action, we need to evaluate the effective potential defined 
in Eq. ([5]) and its derivatives because we use the equations of motion, Eqs. ([7])- ([9]). Because 
there are only two-body interactions in U1 + U2 + U3, this is easily evaluated analytically. On 
the other hand, for path sampling using the conventional transition path sampling (TPS), we 



integrate the overdamped Langevin equation ([3]) for this model from many different initial 
paths. 

Although this is a coarse-grained model, in this paper, we call this model as an "atomistic" 
model or MM DoF because our purpose here is to illustrate our method. The CG DoF 
coupled to this MM DoF will be introduced just below. 

D. CG model for the model polymer and the CG variable 

For the MSES-OM method, we need to introduce CG variables and its OM action. For 
the model polymer introduced above, we take the conventional end-to-end distance z as a 
CG variable, and the "potential" function is introduced as 



t/^^ = a 



f/CG + f;CG //f/CG_^CG 



+ e2 



(18) 



t/CG ^ K^,_,^f^^^^ (19) 



2 

1, 
-/ 
2 



f^2^^ = \HZ-Z2f + V2, (20) 



which is a conventional empirical valence-bond type potential 26|]. In the previous study, it 



is shown that this "single" order parameter can nicely describe the unlooping process of the 



model polymer 30|, so we also use this single variable as the CG DoF. The values of the 



CG parameters are summarized in Appendix B. 

III. RESULTS 

A. Numerical stability of OM dynamics 

It is known that the model polymer described in Sec. Ill CI has two free energy minima with 
respect to the end-to-end distance (this is a basic observation when we build a CG model 
in Sec. Ill Dp . one at the looped configuration, z = 24 nm, and the other at the unlooped 
configuration, 2; = 45 nm. There is a barrier between them, which are characterized by the 
barrier height with ~ 0.4 kcal/mol (looped to barrier top) and ^0.9 kcal/mol (barrier top 
to unlooped). In order to study the unlooping process, the starting point of the OM path 
is fixed at 2; = 24 nm while the end point is made free to move around 2; = 45 nm (see 

10 



Fig. H]). From an independent molecular dynamics run, it has been confirmed that we need 
the trajectory of about 30 fis to get the information on the barrier crossing. Therefore, the 
length of OM path has been taken to be Ny^^^^AtQ^ ~ 40 /is where the OM step size and the 
number of beads are taken as A^om = 10 tg — 16 ns (with t^ being the unit of time scale, 
see Appendix A) and A^bead ~ 2400, respectively. The value of A^q,^ may be taken as small 
as possible, but due to our experience Atgj^ is reasonably small as long as the fluctuation 
of the action is balanced between the bead spring and potential contributions, i.e. each in 
the second term in the right hand side of Eq. (jll). We find the sampling procedure unstable 
when taking AtQj^ too large, say A^qj^ = 1000 tg, when the fluctuation of the potential 
contribution is much larger than that of the bead spring contribution. 

B. Interaction between MM and CG DoF 

Next we show the effect of the coupling between CG and MM paths. In Fig. [21 several CG 
(red) and MM (green) paths during the OM dynamics are drawn with two different coupling 
parameters: A = 10~^ (left) and A = 10~^ (right). We can see that the MM path follows the 
CG path in the strong coupling case (left), while the MM path is more independent from 
the CG path in the weak coupling case (right). 

For a closer analysis on the correlation between MM and CG paths, we define the root 
mean square displacement (RMSD) between the MM and CG DoF as 



5mmcg 



jVCG 
-1 bead / O C 






The histograms for this quantity over 10^ OM steps are shown in Fig. |3]for different coupling 
strengths (A's). We take 6 different values of Aj between 10~^ and 10~^: Aj = 0.01 x 0.25*"^ 
(i = 1, ■ ■ ■ , 6) for 6 different replicas. For the largest A = Ai (=10^^), (5mmcg fluctuates 
around ~ 7, but for the smallest A = Ae (=10~^), the distribution spreads from 5 to 50, 
indicating that the correlation is small between the MM and CG DoF. From this flgure, we 
assume that the paths with Ag is close to those with null coupling (A = 0), because the two 
distributions with A5 and Ae are almost the same. 
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FIG. 2: The end-to-end distance as a function of time (or number of beads) calculated from the 
OM dynamics for MM (green) and CG (red) DoF with the strong (left) and weak (right) coupling 
cases. We overlay 10 "snapshots" from each OM dynamics. There are 10 times difference in the 
MM and CG beads numbers, so every one of 10 beads is shown for the MM DoF. 



C. First passage time distribution and comparison to TPS 



We carried out a MSES-OM calculation with 10^ OM steps. The path data were saved 
every 100 OM steps, and 10'^ data sets were collected. In the MSES-OM calculation, the 
replica exchange was attempted every 10 OM steps, and the acceptance ratio was found 
to be about 30 %. In Fig. HI we show the overlay of snapshots of 100 OM paths from 
the calculation. The time course is indicated by arrows, and the "time difference" between 
the neighboring snapshots is about 4 /is, corresponding to 240 beads. We can see that the 
distribution of paths grows during the time course, but it is hard to tell what is happening 
here. 

In Fig. [5l we show the time evolution of the distribution with respect to the distances 
between "beads" i and j denoted by x^j. For comparison, we also showed the result using the 
conventional OM calculation without Hamiltonian replica exchange, which is called "single 
OM" in the figure. We can see that during the transition the end-to-end distance, x^^, is 
broadly distributed among the two basins (middle-bottom panel). In contrast, the distances 
distributions for x^^ and X35 are already close to those at the final time, indicating that these 
variables, x^g and x^^, first equilibrate before the end-to-end distance does. The difference 
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FIG. 3: The distribution function of the root mean square displacement between the MM and CG 
DoF as defined in Eq. ()2ip for different coupling strengths. Note that we never used Hamiltonian 
replica exchange here. 

between the single OM and MSES-OM methods is found in the distributions for x^^, and X25 
during the transition. This clearly shows the efficiency of our MSES-OM method over the 
single OM method for path sampling. 



Finally we show the first passage time distribution (FPTD) in Fig. [6] 35|, l36|]. This is 
defined as a histogram of the elapsed times in moving from one place to the other. In this 
work we monitor the end-to-end distance z = X15 and measure the time interval where the 
system evolves from z-^^ = 28 to Zg = 35, which correspond to two locations before and after 
the barrier crossing (not free energy minima) [36|. We can see that the barrier crossing 
occurs most likely in 2 /xs, and completes within 35 /is. In order to have a converged result 
for Fig. El we repeated the MSES-OM calculation 10 times from different initial conditions, 
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^-"H 



FIG. 4: Unbiased path ensemble calculated by the MSES-OM method. A hundred of OM trajec- 
tories are superimposed, indicating the fluctuation of the path ensemble. The "time difference" 
between the neighboring snapshots is about 4 /xs, corresponding to 240 beads. 



so the number of path in the ensemble is 10^. 

As a reference, we have also carried out TPS calculations using the Crooks-Chandler 
algorithm 37J for the same system. As can be seen, the FPTD is in reasonable agreement 
with the MSES-OM result. To obtain Fig. El we used 100 independent TPS runs with 10'^ 
TPS steps each. The data is saved every 10 TPS steps, and 10^ data sets are collected 
in total. Thus, the amount of data for TPS calculation and the MSES-OM calculation is 
comparable. The statistical accuracy in Fig. [6] seems to be similar between the MSES-OM 
and TPS methods, at least for the present system. 

However, the comparison of computational effort is difficult in general, since the OM 
methods requires the higher-order derivatives of the potential whereas the parallelized cal- 
culation can be done efficiently. We can at least say that OM method would be of use in 
studying slow and diffusive processes such as conformational changes of large molecules and 
protein folding. There should be much room for further improvement of the OM method 
such as the dominant pathway |20|, the discretization in length [18[, and the use of the 
efficient hessian calculation 381. 
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FIG. 5: Distribution functions for different distances, 3:13, X35 for the upper panel and a;i5,X25 for 
the lower panel, of the model polymer calculated by the single OM and the MSES-OM method. The 
left, middle, and right panels correspond to the distributions at f = 0, 1233AioM, and 2323AtoM) 
respectively. 

IV. CONCLUDING REMARKS 

To efficiently sample huge path space of complex systems, in this paper we proposed a 
multiscale path sampling method based on the Onsager-Machlup (OM) action combined 
with multiscale enhanced sampling (MSES). We applied this MSES-OM method to a model 
polymer with the Asakura-Oosawa interaction, and showed the effectiveness of our method, 
which is comparable to conventional transition path sampling (TPS) with many initial paths. 

We have shown that the MSES method can be easily combined with the OM action. This 
is because the ensemble of stochastic paths subject to overdamped Langevin equation is iso- 
morphic to the canonical ensemble of a polymer of beads connected by harmonic springs. 
Note that the path obtained by the OM or MSES-OM method contains the nonequilib- 
rium information such as the ffist passage time for barrier crossing. This is conceptually 
different from the thermodynamic analysis of barrier crossing, i.e., transition state theory 
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FIG. 6: First passage time distributions (FPTDs) for the model polymer calculated by TPS (green) 
and by the MSES-OM method (red). 

based on the free energy landscape. In order to apply this method to complex molecular 
systems, which is obviously the next step, it is indispensable to make a great use of parallel 
calculation techniques. For a future work, the authors are interested in studying the con- 



formational change of adenylate kinase 



, for which the minimum free energy path was 



recently elucidated by the on-the-fly string method 40 1. 

Another direction to pursue is the Baysian-type inference of the model parameters using 
;he OM action 



31 



41| . It is known that the formulation of the dynamic data assimilation 



42| is almost the same as our MSES-OM method, where the CG DoF corresponds the direct 



variables to be observed, and the MM DoF to the hidden variables of the system studied. 
Hence the MSES-OM formalism should be used for the parameter estimation of general 
dynamical system problems. 
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Finally we mention a possible extension of this approach to quantum systems. Since the 



centroid variable of a quantum particle has been used in TPS 43| and in the string method 



44l |. we expect that we can use the centroid variable in the OM-action based methods as 



well. This will be an extension of the OM path sampling based on the Feynman-Hibbs 



formula which was recently suggested by Faccioli and coworkers 45|, |46 1 . 
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Appendix A: Units for dimensionless variables 

We here summarize the dimensionless units which were used in this paper. We convert 
all the variables (mass, length, energy, time, friction coefficient) into dimensionless ones by 
choosing 

4 4 Q 

mo = p- -ttR^ = -77 X 1.35 g/cm^ x (12.5 nm)^ = 1.10 x 10"^° kg, (Al) 

/o = 10"^ m = Inm, (A2) 

eo = keT = 1.38 x 10"^^ J/K x 300 K, ~ 4.14 x 10'^^ J, (A3) 

to = loj— = 1.63 X 10^9 s = 1.63 ns, (A4) 



eo 
7o = mo/to = 6.75 x lO^^^kg/s. (A5) 

The order of the friction coefficient is given by 

7i = GirriRi = Qrc x 0.05 x 10"^g/cm/s x 12.5 nm 

~ 1.18 X 10"^kg/s = 17570. (A6) 
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If we take as a time step for the Langevin dynamics: At = to/100 = 16.3 ps, which is a 
reasonable choice, the jump magnitude due to the random force is 



^— At = AL ~ 0.01 nm A7 

which seems to be reasonable as well for this model polymer system. 

Appendix B: Parameters for models and simulations 

1. Parameters for the MM DoF 

The parameters used for the model polymer defined in Eqs. ( TT^ -( iT6ll are a = 4nm~^, ei = 
0.24eo, 62 = 70enj_r = 2.5nm, = 0.15, /2j = 12.5 nm, 77 = 5cP,p = 1.35 g/cm , which are 



the same as in [30|] except the value of a. 



2. Parameters for the CG DoF 

The parameters used for the CG model defined in Eqs. (IT8l) - (!20l) are ki = 0.8, /c2 = 
0.02, zi = 24,^2 = 46,t;i = 3.5,^72 = -2.0, e = 5.0, a = 0.7. The parameters for the OM 
action are Atg^ = 100, C^^ = 333, D^^ = 0.003, N^^ = 240 such that the total duration 
becomes the same as the MM case and ksT = 1.0 = D(. (However, this matching of 
the timescales for both DoF might not be necessary.) The value of the diffusion constant 
D = 0.003 was a little bit smaller than that of Fig. 4 in |31)]. In any case, we can estimate 
this value from a short-time simulation. 

Appendix C: On the use of different temperatures for MM and CG DoFs 

At first sight, if we use different temperatures coupled to different parts of the system, the 
system is in a nonequilibrium state which cannot be described by the canonical distribution. 
This is the case in our path ensemble as well as the configuration ensemble in the original 
paper of MSES [26]. However, this problem can be bypassed in the following way. 

Because the formalism is the same for both configuration and path sampling, we em- 
ploy the original setting for configuration sampling. For simplicity we consider only a two- 
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dimensional case with CG variables {Q,P) and MM variables (g,p), but this can be easily 
generalized. 

We set the total (combined) Hamiltonian as 

H = HMM{q,p) + HcG{Q,P) + VMMCG{q,Q) 

= ^ + VMM{q) + ^ + Vcg(Q) + VuMcciq, Q) (CI) 

where m, M are masses for MM and CG DoF. To realize the canonical distribution for 
the total system, we consider the underdamped Langevin dynamics. The equations for the 
original variables are 

-77 = -ImmP ^ A — h \/2'yMMmkBTMMVMM (t) , (C2) 

f = -7ccP - ^ - A^^ + v/27ccMA;,TMMr/cc(t), (C3) 

where 7* is the friction coefficient and ?7*(t) is the Gaussian white noise. Note that we 
used the same temperature Tmm to realize the canonical distribution for the total system 

0^ Q — H/kgTMM 

Changing the CG variables from {Q, P) to {Q, P), which are defined by 

Q = ,/^Q,P = P/,/^, (C4) 

with a = Tmm/Tcg (< l), the Langevin equation for the CG variable becomes 

-TT = -IcqP -^ - A — -~ — + ^^2-fcGMkBTcGVcG{t). (C5) 

at dQ oQ 

Note that the temperature is Tgg in this equation instead of Tmm- 

To have the same a scaling for the terms proportional to A (the third terms in Eqs. (]C2p . 

(]C3p . (ICSP on the right hand side), the interaction Hamiltonian has to be 

^MMCG(g, Q) = \{e{q) - V^Qf (C6) 

where 9{q) represents the transformation from the MM variable to a collective one. 
Hence the above total Hamiltonian can be recast into 

H = HMMiq,p) + ^ + VcGiQ/V^) + V^MMCG(g, Q/V^) (C7) 

with M = M/a. As a result, by simulating Eqs. flC2p and fICSp with different temperatures, 
Tmm and Tcg, we can achieve the canonical distribution for the total system with a single 
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temperature Tmm- In actual calculations, we use only Q space to construct the CG model, 
Vcg(Q) = VcG{Q/\/a), and the interaction is VmhcgI?, Q) = VkMCG(g, Q/V") = f (^(?) - 
QYi it is reasonable to take Vcg(Q) which can mimic the original free energy landscape 
calculated from Vmm{.(1)- 
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